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Abstract 

Tumor heterogeneity is widely considered to be a determinant factor In tumor progression and in particular in its recurrence 
after therapy. Unfortunately, current medical techniques are unable to deduce clinically relevant information about tumor 
heterogeneity by means of non-invasive methods. As a consequence, when radiotherapy is used as a treatment of choice, 
radiation dosimetries are prescribed under the assumption that the malignancy targeted is of a homogeneous nature. In 
this work we discuss the effects of different radiation dose distributions on heterogeneous tumors by means of an 
individual cell-based model. To that end, a case is considered where two tumor cell phenotypes are present, which we 
assume to strongly differ in their respective cell cycle duration and radiosensitivity properties. We show herein that, as a 
result of such differences, the spatial distribution of the corresponding phenotypes, whence the resulting tumor 
heterogeneity can be predicted as growth proceeds. In particular, we show that if we start from a situation where a majority 
of ordinary cancer cells (CCs) and a minority of cancer stem cells (CSCs) are randomly distributed, and we assume that the 
length of CSC cycle is significantly longer than that of CCs, then CSCs become concentrated at an inner region as tumor 
grows. As a consequence we obtain that if CSCs are assumed to be more resistant to radiation than CCs, heterogeneous 
dosimetries can be selected to enhance tumor control by boosting radiation in the region occupied by the more 
radioresistant tumor cell phenotype. It is also shown that, when compared with homogeneous dose distributions as those 
being currently delivered in clinical practice, such heterogeneous radiation dosimetries fare always better than their 
homogeneous counterparts. Finally, limitations to our assumptions and their resulting clinical implications will be discussed. 
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Introduction 

Radiotherapy, the use of ionizing radiation to eliminate 
pathological tissues, is a treatment of choice for more than 50% 
of cancer patients diagnosed with solid tumors [1]. Technical and 
methodological advances have allowed radiation oncology to 
achieve local tumor control in a considerable number of patients. 
However, locoregional recurrence (LRR) remains a problem in 
many clinical settings. For example, a recent study in patients with 
Stage III lung cancer found a 5-year LRR rate of 3 1 % [2] . In 
Glioblastoma Multiforme (GBM), the most common and aggres- 
sive malignant primary brain tumor, LRR approaches 90% [3]. In 
such critical cases, radiotherapy usually results in an initial 
shrinkage of malignancies, followed by a subsequent growth 
recovery that cannot be checked even by resorting to larger 
radiation doses. 

The onset of radioresistance, and its resulting poor prognosis, is 
strongly correlated with the development of significant intratu- 
moral heterogeneity. For that reason, there is growing interest in 
the clinical significance of tumor heterogeneity. In different works 
have been recently demonstrated extensive genetic variations in 



tumor cells due to intratumoral evolution [4], [5]. Moreover, 
tissue-level heterogeneity due to variations in vascular density and 
blood flow has been long since evident in clinical medical imaging. 
In recent years, accumulating evidence suggests that tumor 
heterogeneity is a key factor in the development of therapeutic 
resistance and therefore in radiation therapy outcomes [6], [7], 
[8] . As a consequence, increasing attention is being paid to "dose 
painting" (or "dose sculpting"), a technique which consists in 
prescribing different radiation dosimetries to different regions 
within a given tumor, so that irradiation be boosted in more 
radioresistant (for instance, hypoxic, quiescent, etc.) regions [9], 
[10]. This strategy, which is in sharp contrast with the stiU 
prevailing homogeneous radiation dose delivery approach recom- 
mended by International Commission on Radiation Units (ICRU) 
reports (50, 1993; 62, 1999 and 83, 2010; see [11], [12] and [13] 
respectively), has been made possible by the availability of high- 
precision clinical particle accelerators, and looks particularly 
promising in those cases where current treatment techniques fail 
to provide sufficient tumor control. 

However, in order for dose painting to show its full power, 
detailed information is needed about the internal structure of the 
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tumor to be irradiated. Such information should ideally be 
provided by medical imaging techniques. These however are not 
yet able to distinguish different radiosensitivity regions except in a 
few cases, commonly related to hypoxia effects. On the other 
hand, even a modest miscalculation in the sizes of different 
radiosensitivity subvolumes has been suggested to produce serious 
consequences in cUnical outcomes [14]. In view of current 
technical limitations, the question thus arises of providing tools 
to a) obtain as much information as possible about tumor 
heterogeneity before a radiation dosimetry plan is prescribed, h) 
simulate the effects of dose painting therapies which taike into 
account whatever heterogeneity data are available, and c) compare 
such simulations with those corresponding to standard homoge- 
neous radiation dose distributions currently delivered in clinical 
practice. 

The work herein reported intends to yield some insight into 
these issues. More precisely, in the sequel a mathematical model 
for heterogeneous tumor growth is formulated, and the effects of 
various radiation dose distributions on it are investigated by means 
of computer simulations. Specifically, we consider a situation 
where two tumor cell phenot^pes, cancer cell (CC) and cancer 
stem cell (CSC), are present at an early stage, when the tumor 
consists of about 10^ cells in total. Concerning CCs and CSCs, we 
have assumed that i) CSCs represent only a small percentage of 
the total number of cells at that stage (say, about 15%), ii) CSCs 
have a significantly longer cell cycle duration than CCs and can 
replicate indefinitely, while CCs can perform only a limited 
number of cell divisions, and Hi) CCs and CSCs show quite 
different resistance to radiation, CSCs being more radioresistant 
than CCs. These biological and radiobiological features have been 
reported in the literature, specifically for Glioblastoma Multiforme 
(GBM), where there is mounting evidence of CSCs presence in 
GBM tumors (cf for instance [8], [15], [16], [17]). Growth of the 
heterogeneous tumor thus resulting is simulated by means of an 
agent-based model in which each cell is individually represented 
[18], [19]. Tumor growth is kept track of until a size of 
approximately lO** cells is attained, which roughly corresponds 
to a spheroid of about 1 cubic millimeter in size, a typical volume 
in multi-cellular spheroids (MCS) in vitro growth. At that stage, 
different (homogeneous and heterogeneous) radiation dose distri- 
butions are simulated using the Linear-Quadratic (LQ) model 
[20], [21], and their effects compared. 

An interesting consequence of i) and ii) above is then shown to 
be that, as tumor grows, most of the CSCs concentrate themselves 
within the tumor core, irrespective of their initial distribution at an 
earlier stage. This fact, which will be described to be inversely 
correlated with cell migration rates wlu-n migration is not inhibited 
by cell-cell adhesion (which is the case, for instance, after cells 
undergo an epithelial-mesenchymal transition (EMT) [22]), is then 
used together with Hi) to simulate the effects of different radiation 
dosimetries to achieve tumor control, or in the case where this 
cannot be obtained, to compare tumor heterogeneity (seen as an 
indicator of malignancy in terms of the proportion of CSCs) before 
and after treatment has been delivered. In this context it will be 
shown that for a given amount of radiation, heterogeneous dose 
distributions, where different radiation doses are delivered at 
different regions of the tumor according to the presence of more 
radioresistant cells there, invariably fare better than homogeneous 
ones when sufficient information about tumor spatial heterogene- 
ity is available. In our case, such information will be shown to 
follow from assumptions i) and ii) above. It should be noticed that 
hypotheses i), ii) and Hi) are amenable to experimental validation, 
at least in vitro. 



Our work can be considered as a preliminary step towards 
analyzing preclinical models where larger tumors (of the order of 
cubic centimeters) should be dealt with, several tumor cell 
phenotypes would simultaneously be present (possibly as a 
consequence of mutations) as tumor expands, and vascular 
networks, immune response, and hypoxic and necrotic effects 
are also taken into account. While the case herein considered is 
still far from that situation, the simplicity of the setting selected 
allows us to stress the consi^qiumcc^s derived from tin; minimal 
number of biological and radiobiological assumptions mack' on the 
tumor cell phenot^rpes involved. This last is particularly relevant in 
view of the scarcity of in vivo biological parameter measures 
available. Scaling results up to larger tumor sizes, as well as 
increasing phenotypic and anatomical complexity appear as 
feasible within the same approach, but only after key biological 
data retrievable by non-invasive probing had been identified, and 
their impart on tumor growth elucidated, an objective toward we 
intend to cxmtributc ^vith this work. 

We conclude this introduction by observing that considerable 
attention is being currentiy paid to mathematical modeling as a 
tool towards designing patient-tailored and adaptive therapies; see 
for instance [22], [23], [24], [25], [26], [27], [28], [29], [30] and 
[31]. In partic ular, radiotherapy modeling and simulations have 
been addressed in [32], [33], [34], [35], [36], [37], [38], [39] and 
[40], as well as in [41], [42] and [43] where GBM cases are 
considered. Mathematical models and computer simulations on 
the impact of the presence of CSCs in tumor therapies have been 
discussed in [32], [44] and more recently in [43], where focus is 
made in a GBM case. It is worth to be stressed, however that in the 
cases previously mentioned, the total number of cells simulated 
(and thus the resulting structural complexity) remained way below 
that of the computer simulations arrived at in our current work. 

Materials and Methods 

Tumor Cell Phenotypes Assumptions 

For definiteness model parameter values corresponding to 
Glioblastoma Multiforme (GBM) cell Unes have been used [45], 
[46], [47]. More precisely, we consider a tumor where two 
different phenotypes coexist at an early stage, when we assume a 
preponderant (approximately ?>b% of the total tumor volume) 
proportion of a tumor c:ell phenotype denoted as CC (cancer cell) 
coexisting with a second tumor cell phenotype CSC (cancer stem 
cell), randomly distributed, that roughly represents 15% of the 
total population at that stage. Both phenotypes CC and CSC are 
supposed to possess markedly different biological and radiobio- 
logical properties. In particular, we assume: 

PI.- The duration of cell cycle for CCs is significantiy shorter 
than that of CSCs. In particular, CCs are assumed to divide every 
26 hours. Then, for tumor cell phenotype CSC three cases are 
considered, corresponding respectively to a CSC cycle duration of 
96 hours (four days), 72 hours (three days) and 48 hours (two days). 
Moreover, CCs are assumed to divide a maximum of 15 times, 
while CSCs are able to rephcate indefinitely. 

Concerning property PI, it is currently assumed that CSCs 
proliferate at a slower pace than ordinary cancer cells (see for 
instance [16], [48], [49], [50], [51], [52], [53], [54], [55] and 
[56]). Actually, as observed in the references previously quoted, 
slow-cychng is to be expected from CSCs since such cells belong to 
tumor phenotypes that are highly resistant to current therapies 
(radiotherapy, chemotherapy or combined) and these are targeted 
at killing cycling cells. On the other hand, recent in vivo 
experiments in a mouse model of GUoblastoma to identify and 
isolate CSCs through genetically engineered mice demonstrate the 
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presence of a small pool of slow-cycling and highly tumorigenic 
cells that retain long-term self-renewal ability [55], [57]. We notice 
that cell cycle durations of 24 h for GBM have been reported 
[43], [45], [58], although considerably different cell cycle 
durations, which in particular include the values herein considered 
for CSCs, have been noticed as well [45], [59]. Concerning the 
assumption on the maximum number of CCs replications (see for 
instance [60]), we have selected the value 15 (cf. [32], [61], [62]), 
but our results continue to hold if this number is slightly changed. 
Actually, an arbitrary increase in CSC cycle duration is always 
compatible with our results, as long as CC cycle duration 
continues to be significantly faster. 

In the course of tumor growth, c-ach of the previous tumor cell 
phenotypes may transiendy enter in a quiescent, non-proliferating 
stage, due to contact inhibition. Moreover, replication of CCs is 
always supposed to be symmetric. On the other hand, CSCs will 
be assumed to sustain either symmetric or asymmetric division, in 
which case one CSC and one CC will result from replication. 
Evidence for asymmetric division for CSCs, has been reported in 
[63], [64], [65]. Since reliable estimates about actual probabilities 
of asymmetric division pa do not seem to be available as yet, 
computer simulations will be performed for different choices of 
that model parameter, namely Pa = 0J5, Pa = 0-50 and />a = 0.25 
(cf for instance [32], [61]). 

A second key assumption is: 

P2.- When irradiated, CSCs are significantly more resistant to 
radiation than CCs. 

As a matter of fact, CSCs have been described as a 
comparatively small subpopulation that is highly radioresistant 
[15], [17], [61], [66]. Radioresistance and surviving cell fractions 
are estimated by means of the standard Linear-Quadratic (LQ) 
model [20], [21]. According to it, the surviving fraction of cells 
after a radiation dose D has been delivered, SF{D), is given by: 

5F(Z)) = e-«°^+''^\ (1) 

where D is usually measured in Grays (Gy) (1 Oy is 1 Joule per 
Kilogram), (jx,fi) are the so-called radiosensitivity parameters, 
which depend of the cell phenotype considered, and (J is a 
parameter introduced, as in [32], to distinguish the different 
radiosensitivities of the proliferating and quiescent states for CCs 
and CSCs. Actually, cells in a quiescent state (in the GO cell cycle 
phase) are known to be more resistant to radiation than their non- 
quiescent counterparts [67]. 

It should be noticed that, when estimating the impact of 
radiation according to the LQ_ model, what matters is the 
particular combination of a and P that appears in (1), which 
provides the surviving cell fractions, rather than the separate 
values of a and P by themselves. For definiteness, we take in the 
sequel a = 0A8Gy-\ P = 0.02 Gy-^ and ^^1 = 1.00 for prolifer- 
ating CCs. These radiosensitivity parameters have been reported 
in [47], where in vitro estimates on surviving cell fractions at 
2.0 Of, SF{2), can be found for different GBM cell lines (see silso 
[68]); similar values for a and P have been recently proposed in 
[43]. In particular, 5'i^(2) = 0.36 for proliferating CCs in our case 
(to be compared to the value SF(2) = 0.44 corresponding to a and 
f! parameters proposed in [43]). For quiescent CCs we take 
(Jji=0.85, so that the corresponding value is SF{2) = 0.42. For 
proliferating CSCs the value ^^2=0.30 (5^(2) = 0.73) is taken, 
whereas for rjuiescent CSCs the value (^^2 = 0.20 (5i^(2) = 0.81) is 
selected. Also, such surviving cell fraction ranges at 2.0 Gy have 



been observed and reported in the literature (cf. [43], [68], [69]). 
We point out that the results to be described in this work continue 
to hold when the values selected for the radiosensitivity parameters 

a and p undergo considerable variations, which in particular 
include the ranges considered in the references quoted above. As a 
matter of fact, once assumptions PI and P2 are made, our model is 
shown to be quite robust with respect to changes in its parameters. 

A Three-dimensional {3D) iVlodel of Stochastic Tumor 
Growth 

Different mathematical models of tumor growth and its 
radiation response have been reported in the literature. For 
instance, tumor growth models and radiation effects with 
continuous and discrete populations have been reviewed in [70], 
[71], [72] (see also [33], [34], [35], [36], [73], [74], [75], [76] and 
[77] for more details). On the other hand, the effects of different 
radiation dosimetries have been considered in [32], [33], [34], also 
in [43], [78] for fractionated radiotherapy and in [79] for a case of 
sterc'otactic radiosurgery. Howe\'("r, little seems to be known 
concerning mathematical modeling and computer simulations on 
the effects of heterogeneous radiation dose distributions on 
heterogeneous tumors, as in the case herein examined. 

The model of tumor growth implemented in this work is as 
follows. Within the growing tumor, both tumor cell phenotypes, 
CC and CSC, will be subject to the same kinetic rules. More 
precisely, following [18], [80], [81], a three-dimensional (3D) 
cellular automata (CA) model for tumor growth is developed, 
where each cell is considered as an individual agent. In particular, 
each cell (whether CC or CSC) occupies a single node in a 3D 
unstructured lattice (a lattice with no rotational or translational 
symmetry [81]) thus avoiding symmetry artifacts. Cell division, 
migration, apoptosis (programmed death) and lysis (removal of 
debris) have been included and are represented by stochastic 
processes. Accordingly, each kinetic rule is characterized by a rate, 
and the governing equation to be solved is a multivariate master 
equation (see equation (2)). Figures 1 and 2 show a sketch of the 
processes included and a scheme describing the possible actions 
that a cell is able to perform in the mathematical model 
respectively (for further details, see Document SI provided). 
Nutrient-limited growth is not atxxmnted for in our model. This 
issue, as well as others, could be included at the expense of 
increasing complexity by adding degrees of freedom in the 
computer simulations, but they are not expected to play a 
significant role in a tumor cell colony of the size considered, which 
may be assumed to be fully oxygenated [80]. At any rate, for 
tumors of the size considered in this work our assumption is not 
unlikely. For instance, MIH3T3 cells form tumors of size larger 
than 1 en? without apparent necrotic core even though micro- 
lesions may be observed [80] . 

As to the rules of the model of tumor growth, proliferation is 
only possible for cells located at a node having at least one free 
neighbor in the lattice. In the case that all neighbor sites are 
occupied, a cell enters in a quiescent state due to contact 
inhibition. However, quiescence is abandoned, and cells return to 
their normal state, as soon as one of the surrounding nodes 
becomes free. Proliferation, apoptosis, migration and lysis are 
modeled as stochastic processes occurring with certain rates. A 
Poisson process has been assumed for each individual stochastic 
process and a master equation for the change of the probability of 
the multi-cellular configuration at time t (denoted by the variable 
Z) in terms of the multi-cellular configuration in another state Z' 
at time t' is then used. It reads as follows: 
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Figure 1 . Cell processes mimicked in the model of tumor growth. Schematic representation of the cell processes considered in the model of 
tumor growth (symmetric and asymmetric division, migration, death by radiation, apoptosis (programmed death) and lysis (removal of debris)). 
Notice that CSCs can perform all these cell processes, while replication of CCs is always supposed to be symmetric. See Document SI for further 
details. 

doi:1 0.1 371 /journal.pone.0089380.g001 



' X ' ' = Y.'-z'^zPiZ'AZ"f)-r^^^,P{Zj\Z",t"), (2) 
at ^, 

where p{Z,t\Z" ,t") denotes the conditional probability of finding 
the multi-ceUular configuration Z at time t given the configuration 
was Z" at time t" , fz^z' being the transition rate fi'om 
configuration Z to Z' . Notice that the master equation (2) is a 
balance equation. Indeed, the first term on the right of (2) is a gain 
term that summarizes all transitions that increase the probability of 
finding the corresponding multi-cellular system in configuration Z. 
On its turn, the second term in the right describes transitions that 
move the system away from Z, and thus represents a loss term. 
Equation (2) can be numerically solved if the initial condition 
p(Z,t = Qi) = &{ZQ) is given, where Zq denotes the initial multi- 
ceUular (in our case, tumor cells) configuration. A configuration is 
determined by the spatial distribution of cells and the state of each 
cell (proliferating, quiescent, etc.). In our model both CCs and 
CSCs are able to migrate with the same rate. Migration is 
mimicked by a hopping process allowing any cell to move from 
one lattice site to a free neighbor lattice site. In case several free 
neighbor lattice sites exist, one of them is randomly chosen. In this 
work we have considered two different migration rates (kmig), a 
comparatively low rate (0.025 h~ ') in the range obtained from the 
cell diffusion constant [18] as outlined in the Document SI 
provided, and a higher rate (1.75 '8) as estimated in vitro in [43] 
for a GBM cell line. These will be respectively referred to as low 
and high migration rates in the sequel. It should be stressed that 
we only consider here the case where the motion of a cell from one 
lattice site to another does not depend on the contact energy 
between neighboring cells, but only on the availability of space. In 
that case, the higher the migration rate, the stronger the cell 
dispersion is (see for instance Movies provided as supporting 
information). If however cell-cell adhesion would be considered, 
migrating cells would tend to fill holes and cavities [81], and 



migration will lead instead to tumor compactification. We assume 
that in our current setting this case is substantially included in very 
low migration cases. 

On the other hand, CCs and CSCs undergo programmed cell 
death (apoptosis) (see for instance [32]). Disposal of cellular debris 
resulting from apoptosis is carried out by a lysis process [82], for 
which a lysis rate A:/yj = 0.035/z^' (about 30 A) has been assumed. 
This is about 10-fold less than phagocytosis (digestion of cellular 
debris by macrophages) observed in vivo in [83], but within the 
range reported for in vitro cultures 0.002 h ' for Hybridoma VO 208 
cell line [84] to 0.07 h ' for Fihrobacter succinogenes [85]. 

The master equation (2) has been numerically solved by means 
of the so-called Gillespie algorithm [86], (also called kinetic 
Monte-Carlo algorithm or Bortz-Kalos-Lebowitz algorithm [87], 
see Document SI for more details). Notice that one advantage of 
using a lattice model is the possibility of extending the same 
formalism at larger scales by permitting more than one cell to 
occupy a single lattice site [80]. In order to simulate the resulting 
biological effect when a radiation dose D is delivered, the surviving 
fraction is computed for each tumor cell phenotype according to 
the LQ^ model (1), taking into account the state of each cell, 
proliferating or quiescent. Surviving cells are randomly selected 
out of the total cell number involved. In the next section we wiU 
show a typical starting point (about 10'^ cells) in the computer 
simulations to be described below, as well as the resulting tumor 
once a size of about 10** cells has been reached under the kinetic 
rules just described, and before the radiotherapy treatment is 
applied. For convenience of the reader, we provide in Table 1 all 
parameter values used in our mathematical model to simulate the 
tumor growth and radiation response. 

Results 

As a reference case, we first consider the effect on a fully 
monoclonal tumor (containing only CCs) of a therapeutic 
irradiation protocol, consisting of 30 sessions of 2.0 Gj, each 
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Figure 2. Scheme showing the possible actions that a cell is able to perform in the model. As long as the population size is below a 
prescribed maximum N , it is first tested whether a cell is dead. If so, it undergoes lysis at a certain rate. Alive cells are classified according to CSCs and 
CCs; CCs die and are subject to lysis with a certain rate once they have performed the maximum number of cell divisions prescribed. CCs not having 
yet reached the maximum number of cell divisions and CSCs can undergo apoptosis. Those cells that do not go through apoptosis can migrate if free 
space is available. If they do not migrate and have sufficiently advanced in the cell cycle, they divide. If those cells have not yet reached the end of G2- 
phase, then they continue to progress in the cell cycle. Cells with no free space available at neighboring sites can only progress in the cell cycle. 
Concerning radiation effects, cells are picked randomly and killed according to the corresponding surviving cell fraction estimate. See Document SI 
for details on the technical implementation of the model algorithm. 
doi:1 0.1 371 /journal.pone.0089380.g002 



being homogeneously delivered on the tumor. According to 
standard radiotherapy scheduling, sessions are distributed into 6 
weeks, each week including five sessions fi-om Monday to Friday 
separated by 24 hours intervals, and with a 72 hours interval from 
Friday to next Monday in the following week (with weekend 
interruptions). The total radiation dose delivered with this 
treatment is thus 60 Gy. This is currently considered a standard 
radiotherapy treatment for most GBM tumors [88], [89], [90]. 
The corresponding process is illustrated in Figxire 3, both for the 
cases of low and high migration, under our current assumption 
that migration is not inhibited by cell-cell adhesion. Figure 3 shows 
that tumor prior to treatment grows from week 1 to approximately 
week 5 for the low migration case, and just in one week in the case 
of high migration, until a size of about 1 0 cells is attained. Notice 
that the growth time of the tumor decreases with the migration 
rate due to the decreasing effect of contact inhibition inside the 
tumor. Then radiation therapy treatment starts, and accordingly 
tumor cell number diminishes during the first week (with small re- 
growth between each daily session and weekend interruptions, as 
represented by the knots in the straight line in the plot, see 
Figure 3). The pattern just described is reproduced until tumor 
eradication is achieved at the end of the radiotherapy treatment in 
these cases. 

A heterogeneous tumor containing the two tumor cell 
phenotypes is now considered. Starting from an initial configura- 
tion where 10'^ cells are present, out of which approximately 85% 
are CCs and 15% are CSCs, tumor growth is allowed until a size 
of about 10^ cells is reached (see Figures 4 and 5). Then the impact 
of homogeneous and heterogeneous radiation dose distributions is 
modeled, and computer simulation results are compared in the 
cases where asymmetric division probabilities pa for CSCs are 
Pa = 0.75, p„ = 0.50 and Pu = 0.25, the CSC cycle duration is taken 
to be 96 A, 72 h and 48 h, and the low and high migration rates are 
considered. The results obtained wlU show that the standard 



irradiation protocol described before fails now to achieve tumor 
control in any of the cases considered. To compare the dynamics 
of the tumor resulting after irradiation with respect to its pre- 
treatment stage, computer simulations are stopped once the pre- 
treatment population size of about lO'' cells is again obtained. 

As shown in Figure 4, for = 0.25, CSC cycle duration equal 
to 96 h and the low migration rate, the more radioresistant tumor 
cell phenotype CSC is confined within an inner, smaller region 
when irradiation is started. Such spatial CSCs distribution is 
neither a priori imposed, nor a consequence of the specific CSC 
cycle duration or the asymmetric division probability considered 
(see Table 2). It is due instead to the difference of the CSC and CC 
cycle durations. Indeed, a robust emerging feature is now 
observed. Namely, due to asymmetric division CSCs produce a 
certain fraction of CCs. Both CCs and CSCs then compete for 
resources including free space at the tumor border [18], [91]. For 
sufficiently small micro-motHity, that competition is controlled by 
cell replication. As CCs proliferate faster than CSCs, they have a 
selective advantage in the competition for free space and wiU 
eventually outcompete the CSCs in the border region of the 
tumor, if (as it happens in our case) to achieve such dominance less 
replications are needed than the maximum number that CCs can 
perform. The precise timing depends on the relation of the cell 
cycle duration for CSCs vs. CCs, Pa, and the fraction of CSCs in 
the initial population at 10'^ cells (notice that this fraction would 
itself be determined by />„ and CSC cycle duration if the 10'^ cells 
would already have emerged by replication from a single initial 
CSC). 

Therefore, for a low migration rate, CSCs wiU be contact- 
inhibited by the fast proliferating CCs. As a consequence, CSCs 
will remain confined in an inner region in that case. Actually, on 
assuming a cell diameter of 20 \lm, the diameter of the tumor in all 
cases is then of about 2680 \lm (with a standard deviation of 56 \m 
over 20 simulations performed for each parameter set considered) 
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Table 1. Model parameters used in computer simulations of tumor growth and radiotherapy treatments. 



Description 


Symbol 


Value/Range 


Source 


Migration rate 


kmig 


0.025 h"Vl.75 h"' 


[18], [43] 


Apoptosis rate 




4.17x10"" h"' 


[32] 


Lysis rate 


^lys 


0.035 


(Assumed) 


Radiosensitivity (LQ model) 


1 


0.48 6y"' 


[47] 


Radiosensitivity (LQ mode!) 


j8 


0.02 Gy"^ 


[47] 


Radiosensitivity (LQ model}: Proliferating (CC) 




1.00 


[68], [69] 


Radiosensitivity (LQ model}: Quiescent (CC) 




0.85 


[68], [69] 


Radiosensitivity (LQ model): Proliferating (CSC) 


ip2 


0.30 


[68], [69] 


Radiosensitivity (LQ model): Quiescent (CSC) 




0.20 


[68], [69] 


CC cycle duration 


■ta 


26 h 


[45] 


CSC cycle duration 




48 h, 72 b, 96 h 


(Assumed) 


Asymmetric division probability (CSC) 


Pc 


0.75, 0.50, 0.25 


(Assumed) 


Maximum number of cycle divisions (CC) 




15 


[32], [61] 



Values for those parameters not found in the literature were assumed (see detailed explanation for lysis rate). In the remaining cases (asymmetric division probability 
and CSC cycle duration) some values were assumed, and the impact of different parameter sets on the resulting effects was subsequently analyzed. 
doi:l 0.1 371 /journal.pone.0089380.t001 




Figure 3. Standard radiotherapy treatment in a homogeneous tumor for the low and high migration cases. Cell growth curves are 
shown corresponding to homogeneous tumor growth for the low and high migration cases when only CCs are present (see respectively (A) and (C)). 
Tumor growth is allowed unchecked from a size of about 10^ cells until about 10* cells are present, which approximately occurs at day 30 
(respectively at day 7) since the beginning of the process. Then, a homogeneous treatment corresponding to 30 sessions of 2.0 Gy each is delivered. 
In all cases, sessions are scheduled along 6 weeks separated by 24 hours intervals except for weekends, where a 72 hours interval is allowed. 
Radiotherapy treatment is thus completed 40 days afterwards its beginning (about 70 and 47 days since the initial stage respectively). Data 
corresponding to 20 simulations (with different seeds of a random number generator) are presented. Notice that the vertical coordinate is 
represented in a logarithmic scale. In (B) and (D) tumor stages are represented when radiation therapy is started (about 10* cells in total) for the low 
and high migration cases respectively. Depicted in dark and light green are proliferating and quiescent CCs. Dead cells are represented in black. See 
Movies SI and S2 for an example of simulations represented in (A) and (C). 
doi:10.1371/journal.pone.0089380.g003 
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Figure 4. Simulated growtli of a heterogeneous tumor with the low migration rate. Depicted in darl< and light green (respectively, dark 
and light red) are proliferating and quiescent CCs (respectively, proliferating and quiescent CSCs). Dead cells are represented in black. (A) An initial 
stage where about 10^ cells, distributed into tumor cell phenotypes CC (85%) and CSC (15%), are present. (B) Tumor stage when radiation therapy is 
started (about 10^ cells in total). In the middle image, the location of the inner region where 100% of CSCs are concentrated is shown for the case 
when Pa = 0.25 and CSC cycle duration equal to 96 h. A 3D transversal cut is performed in the middle of solid figures (A) and (B) (left), so that its 
interior could be seen (middle and right) respectively. (C) Representation of the transversal cut showed in (B) for a slice of two cell diameters. Notice 
the little space existing between cells. 
doi:1 0.1 371 /journal.pone.0089380.g004 



and the volume of this inner region where 100% of CSCs are 
located varies from the 15% to 25% of the total tumor volume, 
when asymmetric division probability and the CSC cycle duration 
are allowed to change (see Table 2, and Figures in the Document 
SI provided). 

On the other hand, when the high migration rate is considered, 
CSCs are not fully concentrated in an inner region of the tumor 
(see Figure 5 for /i^ = 0.25 and CSC cycle duration equal to 48 h). 
However, we can define an inner region where at least 80% of 
CSCs are located (see Figure 6). When asymmetric division 
probability and CSC cycle duration are allowed to change in the 
parameter range considered, this inner region approximately 
represents between 21% to 40% of the volume where 90% of cells, 
both CCs and CSCs, are located (see Table 3, and Figures in the 
Document SI provided). In this case, the diameter of the tumor for 
all cases is about 5294 [im (with a standard deviation of 778 |a,m 
over 20 simulations performed for each parameter set considered). 
In Tables 2 and 3 the number of CSCs just before treatment starts 
is shown, so that its dependence with migration rate, asymmetric 
division probability and CSC cycle duration can be observed. 
Actually, the number of CSCs existing before treatment starts is a 



key factor to estimate tumor resistance to radiation therapy, as we 
win recall below. 

Bearing these facts in mind, it turns out that tumor control can 
be obtained in all cases when a radiation boost is applied at such 
internal regions. More precisely, in the case of low migration we 
observe that tumor control can be achieved for CSC cycle 
durations equal to 96 h and 72 h, when 2.5 Gy (for the case 
Pa = 0.75), 2.9 Gj (for p„ = 0.50) and 3.3 Gy (for ^,, = 0.25) are 
respectively delivered within the largest inner sphere containing 
100% of CSCs, and 2.0 Gy is delivered in the rest of the tumor, 
according to the former standard fractionation protocol (5 days a 
week along 30 sessions at 24 hours intervals except for weekends). 
However, when CSC cycle duration is 48 h, tumor eradication is 
not possible with these heterogeneous therapies under the same 
conditions. In that case, to obtain tumor control, the dose 
delivered in the inner region has to be raised to 2.7 Gji, 3.4 Gy and 
3.9 Gy respectively (see Table 4 and the Document SI provided). 
Notice that these radiation doses increase as asymmetric division 
probability decreases. 

Let us now consider the same heterogeneous therapies for the 
case of high migration. We now select an inner region where 80% 
of CSCs are located. Considering these heterogeneous therapies 
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Figure 5. Simulated growtli of a heterogeneous tumor with the high migration rate. Depicted in darl< and light green (respectively, dark 
and light red) are proliferating and quiescent CCs (respectively, proliferating and quiescent CSCs). Dead cells are represented in black. (A) An initial 
stage where about 10^ cells, distributed into tumor cell phenotypes CC (85%) and CSC (15%), are present. (B) Tumor stage when radiation therapy is 
started (about 10* cells in total). In the right image, the spatial distribution of CCs and CSCs is shown for the case when p„ = 0.25 and CSC cycle 
duration equal to 48 h. A 3D transversal cut is performed in the middle of solid figure (B) (left), so that its interior could be seen (right). (C) 
Representation of the transversal cut showed in (B) for a slice of two cell diameters. Notice the comparatively large (with respect to Figure 4) space 
observed between cells. 
doi:1 0.1 371 /journal.pone.0089380.g005 



for each case as before, tumor control is now obtained only for 
= 0.75 with CSC cycle durations of 96 h and 72 h. This is due 
to the fact that i) the high migration rate permits less contact 
inhibition, which in turn allows for rapid re-growth, and ii) there 
are about 20% of CSCs which only receive a radiation dose of 
2.0 Gy. Therefore, to obtain tumor control it is not only necessary 
to increase tiie radiation dose in the inner region, but also in the 
outer one (see Table 5 and the Document SI provided). The 
radiation doses of the heterogeneous therapies required to obtain 
tumor control are provided for each case of migration rate 
considered in Tables 4 and 5. The respective temporal evolution of 
the number of each tumor cell phenotype is shown in Figure 7 (A), 



(C), (E) and (G) for different values of asymmetric division 
probability, migration rate and CSC cycle duration. 

It may appear at first glance that the successful results obtained 
for heterogeneous dosimetries could be a consequence of the 
overall radiation dose delivered over the tumor being larger than 
that administered according the standard irradiation protocol 
(2.0 Oy a day, 5 days a week at 24 hours intervals, with weekend 
interruptions and 60 Gy in total). However, heterogeneous 
dosimetry turns out to be crucial to achieve tumor control. In 
particular, tumor control fails to be attained when we deliver an 
averaged homogeneous dose (AD), corresponding to the same 
global radiation energy as in the heterogeneous dosimetry, carried 
out along a similar scheduling. The corresponding AD is given by: 
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Table 2. Estimates of the tumor inner region diameter and 



number of CSCs before irradiation for the low migration case. 






Pa =0.75 




Pa =0.50 




Pa = 0.25 






Diameter 


CSCs 


Diameter 


CSCs 


Diameter 


CSCs 


96 h 


1426.3 /im 


16871 


1488.6 /im 


18718 


1513.5 /im 


20448 




[41.07 /im] 


[56.12] 


[43.24 /im] 


[71.33] 


[53.70 ;im] 


[77.30] 


72 h 


1473.6 /im 


17125 


1539.3 /im 


19366 


1587.8 /im 


21092 




[40.85 /im] 


[67.50] 


[58.67 /im] 


[102.07] 


[40.45 ;im] 


[85.48] 


48 h 


1498.1 /im 


17785 


1595.9 /im 


19829 


1682.2 /m 


21953 




[47.21 /im] 


[62.45] 


[73.16 ;im] 


[88.98] 


[47.74 ;im] 


[140.16] 





Diameter is that of an inner sphere where 100% of CSCs are located. CSCs 
number is computed before radiation therapy treatment starts. Within brackets 
the corresponding standard deviations are also provided. Data corresponding 
to 20 simulations (with different seeds of a random number generator) for each 
case considered. See also Figures in the Document SI for further details. 
doi:1 0.1 371 /journal.pone.0089380.t002 



where Vont are the volume of the internal sphere and the 
remaining shell considered; Do^t are the radiation doses 

delivered over the internal and external regions just described, and 
Vtot is the total volume of the tumor. In the case of low migration 
the inner region is that where 100% of CSCs are located (see 
Table 2 for values of the diameter of this inner region for each 
case) and the volume of the outer region is computed with respect 
to the average diameter of the tumor at the beginning of the 
radiotherapy treatment (2680 \urt). However, for the case of high 
migration the inner region is now selected as that where 80% of 
CSCs are located (see Table 3 for further details). Indeed, since 
some cells are now isolated far from the tumor bulk, instead of 
defining the tumor radius as the distance from its center of mass to 
the farthest cell, to compute the averaged dose the volume of the 
tumor is now considered as that of the region where 90% of total 
cells (CCs and CSCs) are located, where the diameter is about 
3 1 20 |Xm (with a standard deviation of 1 86 m over 20 simulations 
performed for each parameter set considered). The reason for this 
assumption is that it will yield a higher AD than that obtained 
when considering the maximum diameter of the tumor (5294 |a,m), 
which will extend to regions sparsely occupied by tumor cells. 
Hence the averaged homogeneous therapies thus derived will 
deliver higher radiation doses than those that would be obtained if 
the outer shell were defined as that where all tumor cells are 
contained. 

The averaged dose [AD) per session according to equation (3) is 
shown in Tables 4 and 5, for each of heterogeneous therapies 
described before. These AD vary from 2.10 Gy to 2.32 Gy for the 
low migration case and from 2.10 Gylo 2.52 Gy in the case of high 
migration for the asymmetric division probabilities and CSC cycle 
durations considered. Notice that the total radiation doses 
delivered by these averaged homogeneous therapies are higher 
than 60 Gy (the value corresponding to the standard irradiation 
protocol) for all cases. The total radiation doses corresponding to 
these new dosimetries range between 63.0 Gy to 69.6 Gy for the 
case of low migration and 63.0 to 75.6 Gy for the high 
migration case. Some of these results are illustrated in Figure 7 (B), 
(D), (F) and (H). On the other hand, in Figure 8, further details of 
the time evolution of the tumor colony are provided during and 



after an homogeneous radiation therapy delivering an 
AD = 1.\^ Gy for p„ = 0.75 and AD = l.l'h Gy for p„ = 0.25 with 
the low migration rate and CSC cycle duration of 96 h. The cases 
Pa = 0.25 with the high migration rate and CSC cycle durations of 
96 h and 48 h, an ^Z) = 2.36 Gy and AD = 252 Gy respectively 
are also included. Notice that in the case of high migration 
(Figwe 8 (B)) more cells remain isolated at the end of the treatment 
compared with the case of low migration (Figure 8 (A)). In that 
Figure, when tumor control is not achieved, computer simulations 
are performed until the sumving tumor reaches a size approx- 
imately equal to 10*" cells, the number of cells it had before 
radiotherapy started. 

Tables 4 and 5 reveal that tumor recurrence occurs in all cases 
for a homogeneous therapy delivering the corresponding average 
dose {AD). Besides, the number of CSCs in the tumor at the end of 
the radiotherapy treatment decreases with Pa and CSC cycle 
duration (see Figure 9, and Tables in the Document SI provided). 
In the case of low migration, for the heterogeneous therapies 
failing to achieve tumor control, the number of CSCs remaining 
alive at the end of the recurrence tumor stage is 107, 1785 and 
4457 respectively, with the corresponding standard deviations 
being 8.53, 78.31 and 232.67 (see Figure 9 (A) to compare with the 
corresponding averaged homogeneous therapies). These values 
correspond to the cases Pa = 0.75, /i^ = 0.50 and pa = 0.25 with a 
CSC cycle duration of 48 h. In Figure 9 (B), the number of CSCs 
at the end of the recurrence tumor stage is provided in the case of 
high migration for the heterogeneous therapies delivering 2.5 Gy 
(for the case p„ = 0J5), 2.9 Gy (for Pa = 0.50) and 3.3 Gy (for 
Pa = 0.25) in the inner region, and 2.0 Gy in the rest of the tumor. 
Notice that, even when tumor control cannot be achieved with the 
heterogeneous therapies, the corresponding averaged homoge- 
neous therapies always have more CSCs at the end of the 
recurrence tumor stage (see Figure 9 (C)). Moreover, in some cases 
that number of CSCs is larger than before the treatment started, 
resulting in more radioresistant tumors after treatment (see 
Document SI for further details). 

Thus, to achieve full eradication of a tumor consisting of two 
different tumor cell phenotypes, heterogeneous dosimetry is 
crucial. Actually, the choice of a minimal radiation dose sufficient 
to achieve tumor control depends on the value ofpa, the CC and 
CSC cycle durations and on the internal spatial distribution of 
CSCs. In Tables 4 and 5, we describe the heterogeneous radiation 
therapies needed to achieve tumor control, and the corresponding 
averaged homogeneous therapies are also provided. Interestingly, 
the corresponding averaged homogeneous therapies in each case 
fail to obtain tumor control (see also Tables in the Document SI 
provided). Moreover, homogeneous therapies needed to obtain 
tumor control are also provided in Tables 4 and 5. One readily 
sees that in all cases higher total radiation doses are needed for 
homogeneous than for heterogeneous therapies. 

On the other hand, considering that for all choices of model 
parameters the AD is higher than 2.0 Gy, this implies that tumor 
recurrence will also occur for the standard irradiation protocol 
(2.0 Gy a day, 5 days a week at 24 hours intervals, with weekend 
interruptions and 6. Gy in total) for each case ofpa, migration rate 
and CSC cycle duration considered. In terms of the number of 
remaining CSCs after treatment is completed, recurrence is 
certainly weaker when AD is delivered than for the standard 
fractionation protocol, as one could expect from the comparative 
increase in radiation delivery. Moreover, tumor control cannot be 
achieved for each case of Pa, migration rate and CSC cycle 
duration considered even when the homogeneous therapy 
delivering the average radiation dose is rescheduled in 7 days a 
week along 30 sessions at 24 hours intervals, without weekend 
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Inner region (gQo/^ of ccs + CSCs) 

(80% of CSCs) 





Figure 6. Spatial distribution of CSCs for a heterogeneous tumor witKi the high migration rate. From left to right tumor stage when 
radiation therapy is started (about 10^ cells in total) with the high migration rate, 3D transversal cut in the middle of the tumor, region where 80% of 
CSCs are located (yellow) and region where 90% of total cells (CCs and CSCs) are located (yellow and blue). (A) For the case /i„ = 0.75 and (B) for 
Pa = 0.25 considering CSC cycle duration equal to 48 h. Depicted in dark and light green (respectively, dark and light red) are proliferating and 
quiescent CCs (respectively, proliferating and quiescent CSCs). Dead cells are represented in black. 
doi:1 0.1 371/journal.pone.0089380.g006 



interruptions (see Figure 10 for some examples of averaged 
homogeneous therapies with this fractionation protocol). 

Since in many clinical scenarios radiation doses are mostly 
limited by damage inflicted at neighboring organs at risk and 
healthy tissues (see [92], [93] and [94]), it is important to estimate 
what amount of tumor control can be achieved when radiation 
dose distributions are kept as low as possible. In what follows, we 
shall consider a heterogeneous therapy for which the average 
radiation dose is approximately equal to 60 Gy, and a case of 
hyperfractionation (a type of scheduling consisting of compara- 
tively many sessions, usually more than 1 per day, with low 
radiation doses [95]); cf [96] for a specific study on GBM tumors. 
We shall see that in these cases a heterogeneous radiation dose 
distribution also yields better results than its averaged homoge- 
neous equivalent, even when tumor control is not achieved. 



Table 3. Estimates of the tunnor inner region diameter and 
number of CSCs before irradiation for the high migration case. 









Pa = 0.75 




Pa =0.50 




Pa =0.25 






Diameter 


CSCs 


Diameter 


CSCs 


Diameter 


CSCs 


96 h 


1857.4 fim 


20454 


1 986.9 fim 


27916 


2035.6 fim 


35087 




[74.72 ^lm] 


[256.53] 


[51.30 fim] 


[811.67] 


[77.80 fim] 


[1066.54] 


72 h 


1906.2 fim 


21178 


2043.1 fim 


28847 


2158.7 fim 


37686 




[54.46 fim] 


[322.63] 


[78.51 fim] 


[861.47] 


[94.21 fim] 


[859.31] 


48 h 


1983.8 fim 


21944 


2139.3 fim 


30119 


2294.8 fim 


41629 




[69.64 fim] 


[506.92] 


[81.17 fim] 


[872.98] 


[62.60 fim] 


[1040.65] 



Diameter is that of an inner sphere where 80% of CSCs are located. CSCs 
number is computed before radiation therapy treatment starts. Within brackets 
the corresponding standard deviations are also provided. Data corresponding 
to 20 simulations (with different seeds of a random number generator) for each 
case considered. See also Figures in the Document SI for further details. 
doi:l 0.1 371/journal.pone.0089380.t003 



Consider first the case of low migration, CSC cycle duration equal to 
96 h and where the value of the total radiation dose is a bit less than 
60 Gy for heterogeneous therapies consisting of 2.3 Gy for pa = 0J5, 
Pa = 0.50 and Pa = 0.25 within tlie largest inner sphere containing 
100% of CSCs and 1.8 Gj; in the rest of the tumor delivered 5 days a 
week along 30 sessions at 24 hours intervals with weekend 
intenoiptions. Computer simulations show that these radiation 
dosimetries fare better than their averaged homogeneous versions, 
where total AD lies between 56 Gy and 58 Gy for tlie values of Pa 
considered (see Figure 1 1 (A), (B) and (C)). On the other hand, similar 
results can be obtained for a lower total radiation dose when the time 
lapse between sessions Ls also shortened. More precisely, consider the 
same cases but now for heterogeneous therapies consisting of 1.7 Gy 
within the largest inner sphere containing 100% of CSCs and 1.2 in 
tlie rest of the tumor, delivered in two sessions per day, 5 days a week 
along 30 sessions 1 2 hours intervals, with weekend interruptions. While 
tumor control is not achieved, tumor radioresistance, measured in 
terms of the final proportion of CSCs, turns out to be lower for 
heterogeneous dosimetries than their averaged versions (see Figure 1 1 
(D), (E) and (F)). Notice that in this case, the total averaged doses 
delivered by the heterogeneous dosimetries considered are much 
smaller than 60 Gy (between 39 Gy and 39 Gc for the values of Pa 
considered). 

Figure 1 1 shows that there is tumor recurrence in aU cases. 
However, it turns out that the number of CSCs remaining at the end of 
the recurrence tumor stage after radiation therapy is lower than that 
existing prior to therapy in all cases. Thus tumors surviving this therapy 
can be considered as less radioresistant than they were before radiation 
therapy started. An inspection of Figure 1 1 (C) and (F) quickly shows 
that in our current cases heterogeneous therapies yield better results 
than its averaged homogeneous counterparts previously discussed. For 
completeness, we provide estimates on the total number of CSCs after 
treatments are concluded and recurrence appears for each case of Pa 
considered (see Fignre 11). 

We conclude this Section by remarking on the dependence of 
our model of tumor growth, and the results derived from its 
analysis, on data and parameters assumed. 
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Table 4. Classification of heterogeneous an 


d homogeneous 


radiation therapies for the low migration case. 






Heterogeneous therapy 




Homogeneous therapy 




Pa Tcsc No Control 


Control 


No Control 


Control 


0.75 96 h 


2.0 Gy-2.5 Gy"' 


2.10 Gy'" 


2.5 Gy 




[63.0 Gy] 


[63.0 Gy] 


[75.0 Gy] 


72 h 


2.0 Gy-2.5 Gy'^' 


2.10 Gy'^' 


2.5 Gy 




[63.0 Gy] 


[63.0 Gy] 


[75.0 Gy] 


48 h 2.0Gj'-2.5Gy 


2.0 Gy-2.7 Gy'"' 


2.10 Gy'^'/ 
2.12 Gy'"' 


2.7 Gy 


[63.0 Gy] 


[63.6 Gy] 


[63.0 Gy]/ 
[63.6 Gy] 


[81.0 Gy] 


0.50 96 h 


2.0 Gy-2.9 Gy'^' 


2.15 Gy'=' 


2.9 Gy 




[64.5 Gy] 


[64.5 Gy] 


[87.0 Gy] 


72 h 


2.0 Gy-2.9 Gy"^' 


2.17 Gy"^' 


2.9 Gy 




[65.1 Gy] 


[65.1 Gy] 


[87.0 Gy] 


48 h 2.0Gj'-2.9G3''" 


2.0 Gy-3.4 Gy"'' 


2.19 Gy'"/ 
2.30 Gy"" 


3.4 Gy 


[65.7 Gy] 


[69.0 Gy] 


[65.7 Gy]/ 
[ 69.0 Gy] 


[1 02 Gy] 


0.25 96 h 


2.0 Gy-3.3 Gy'" 


2.23 Gy'" 


3.3 Gy 




[66.9 Gy] 


[66.9 Gy] 


[99.0 Gy] 


72 h 


2.0 Gy-3.3 Gy"°' 


2.27 Gydo) 


3.3 Gy 




[68.1 Gy] 


[68.1 Gy] 


[99.0 Gy] 


48 h 2.0G;'-3.3G/"' 


2.0 Gy-3.9 Gy"^' 


2.32 Gy'"'/ 
2.47 Gy"^' 


3.9 Gy 


[69.6 Gy] 


[74.1 Gy] 


[69.6 Gy]/ 
[74.1 Gy] 


[117 Gy] 



In all cases, treatment sessions were scheduled along 6 weeks separated by 24 hours intervals except for weekends, where a 72 hours interval is allowed. Data 
corresponding to 20 simulations (with different seeds of a random number generator) are presented. In the heterogeneous therapies, radiation doses are specified both 
for the outer (left) and inner (right) tumor regions, each case being indexed from (1 ) to (12). The averaged dose for any of the previous cases is labeled with the same 
number in the columns corresponding to homogeneous therapies. Within brackets the total dose of the radiation therapy treatment is also provided. See Tables and 
Figures in the Document SI for further details. 
doi:1 0.1 371 /journal.pone.0089380.t004 



To begin with, our results are not restricted to the figure 
selected (15%) for the proportion of GSCs within the tumor at the 
initial stage. In fact, they continue to hold as long as the more 
radioresistant tumor cell phenotype CSC represents a small 
percentage of the total tumor cell count. A particularly interesting 
hmit case is that when tumor growth starts from a single CSC. 
Then for each value of^a (0.75, 0.50, 0.25), CSC cycle duration 
(96 A, 72 h and 48 A) and migration rate k^ig (0.025 A"', 1.75 /;"') 
considered, the number of CSCs present when tumor has reached 
a size of about 10 cells (just before radiation treatment starts) is 
much smaller than that corresponding to the cases considered in 
this work. For instance, in the case of CSC cycle duration equal to 
48 h. Pa = 0.25, for the low and high migration cases the number 
of CSCs before the treatment starts is about 5956 and 14316 (with 
standard deviations of 129 and 530 over 20 simulations performed 
for each parameter set considered) respectively. The correspond- 
ing values for the case considered in this work are 21953 and 
41629 respectively (see Tables 2 and 3). Moreover, the internal 
region where CSCs remain confined is smaller (1120 Urn and 
1840 |Xm with standard deviations of 39.7 [Im and 71.4 [im over 20 
simulations performed for each parameter set considered, respec- 
tively) than that reported in Tables 2 and 3. Hence any 



radiotherapy treatment that achieves tumor control in our case 
also does so for tumors staring from a single CSC under 
assumptions above. 

On the other hand, we have made use of the assumption that 
the duration of cell cycle for CSCs is significandy longer than that 
of CCs, a hypothesis commonly assumed in the literature (cf [16], 
[48], [49], [50], [51], [52], [53], [54], [55] and [56]). This fact 
notwithstanding, our model can be used to examine also the 
opposite situation, that is the case where CSC cycle duration is 
equal or smaller than that of ordinary CCs. As an example, we 
have considered the cases where CSC cycle lasts 26 hours 
(respectively 18 hours), which is equal to (respectively less than) the 
26 hours cell cycle selected for CCs. As one can expect, the inner 
core where most CSCs remain concentrated is now larger than 
when slow-cycling CSCs is assumed. In particular, in the case of 
low migration and for a CSC cycle duration of 26 h, such internal 
regions (where 100% of CSCs are located) range from 20% to 
83% of the total tumor volume for values of^a equal to 0.75, 0.50 
and 0.25. Besides, when CSC cycle duration is taken to be 18 h 
that internal volume further expands, ranging now between 23% 
and 100% of the total tumor volume. Additional details, including 
the number of CSCs present when tumor size reaches about 10^ 
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Table 5. Classification of heterogeneous an 


d homogeneous radiation therapies for the high migration case. 








Heterogeneous therapy 




Homogeneous therapy 




Pa ■''esc 


No Control 


Control 


No Control 


Control 


0.75 96 h 


- 


2.0 Gy-2.5 Gy'" 


2.5 Gy<" 


2.5 Gy 






[63.0 Gy] 


[63.0 Gy] 


[75.0 Gy] 


72 h 




2.0 Gy-2.5 Gy'^' 


2.11 Gy'^' 


2.5 Gy 






[63.3 Gy] 


[63.3 Gy] 


[75.0 Gy] 


48 h 


2.0 Gy-2.5 Gy'^' 


1 -I -I -T ,(4) 

2.2 Gy-2.7 Gy 


2.1 3 Gy'^'/ 
2.33 Gy"" 


2.7 Gy 
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In all cases, treatment sessions were scheduled along 6 weeks separated by 24 hours intervals except for weekends, where a 72 hours interval is allowed. Data 
corresponding to 20 simulations (with different seeds of a random number generator) are presented. In the heterogeneous therapies, radiation doses are specified both 
for the outer (left) and inner (right) tumor regions, each case being indexed from (1) to (16). The averaged dose for any of the previous cases is labeled with the same 
number in the columns corresponding to homogeneous therapies. Within brackets the total dose of the radiation therapy treatment is also provided. See Tables and 
Figures in the Document SI for further details. 
doi:1 0.1 371 /journal.pone.0089380.t005 



cells and the case of high migration, are provided (see Document 
SI). 

A case which has not been addressed in this work is cancer cell 
plasticity, a hypothesis that has been advanced to better 
understand the onset of resistance after therapy; see for instance 
[97], [98] and [99]. According to this scenario, in addition to 
CSCs giving raise to CCs by asymmetric division, a (supposedly 
small) percentage of CCs may transform to a CSC phenotype, 
possibly as a reaction to radiation therapy. Although little 
quantitative information about this process seems to be available 
as yet, including such type of process in our model is possible. To 
support this statement, we have studied a particular example. 
Specifically, we have examined a model situation where a small 
percentage of CCs are transformed to CSCs along the radiation 
treatment. As expected, any increase in the number of CSCs 
results in increased malignancy, measured in terms of higher 
resistance to radiation therapy. However, our conclusion that 
heterogeneous, tumor-adapted radiation therapies fare better than 
their corresponding averaged homogeneous versions continues to 



hold. Details on this study can be found in the Document SI 
provided. 

Discussion 

Tumor heterogeneity is being increasingly recognized as a key 
obstacle to achieve successful tumor control, either by means of 
radiotherapy, chemotherapy or through the use of combined 
therapies. Indeed, it is well known that tumors at an advanced 
stage contain different tumor subpopulations, which might have 
been generated as a consequence of sequential mutations of one 
initial clonogenic line, or could result from the presence of Cancer 
Stem Cells. Moreover, it is expected that such cell phenotypes may 
considerably dififer in their biological and radiobiological proper- 
ties, and in particular in their resistance to radiation (cf for 
instance [1,5], [6], [60], [8], [100], [101]). 

Accordingly, it has been proposed that the clinical prognosis of a 
given tumor would critically depend on the information that may 
be gathered about its internal heterogeneity, and more precisely. 
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Figure 7. Comparing heterogeneous and averaged homogeneous radiation therapies in a heterogeneous tumor for different 
model parameters. Cell survival curves for 20 simulations (with different seeds of a random number generator) in the cases = 0.75 and = 0.25 
for CSC cycle durations equal to 96 h and 48 h with the high and low migration rates are shown. The time evolution for CCs and CSCs is represented 
in green and red respectively. (A, C, E, G) Results for heterogeneous therapies consisting of 2.5 Gy and 3.3 Gy in the inner sphere and 2.0 Gy in the rest 
of the tumor. (B, D, F, H) Results for the related averaged homogeneous therapies corresponding to 2.10 Gy, 2.23 Gy, 2.36 Gy and 2.52 Gy 
respectively. (A, B, C, D) Results for the cases p„ = 0J5 and /)„ = 0.25 with the low migration rate and CSC cycle duration equal to 96 h. (E, F, G, H) 
Results for the case /7„ = 0.25 with the high migration rate and CSC cycle durations equal to 93 h (E, F) and 48 h (G, H). In all cases 30 sessions are 
scheduled along 6 weeks, separated by 24 hours intervals except for weekends, where a 72 hours interval is allowed. Radiation is applied when the 
total cell count is about 1 0* cells. Notice that the vertical coordinate is represented in a logarithmic scale. See Movies S3, S4, S7 and S8 for an example 
of simulations represented in (B), (D), (G) and (H) respectively. 
doi:1 0.1 371/journal.pone.0089380.g007 
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about the identification of regions within it with different 
sensitivities to a given therapy. In principle, once the spatial 
distribution of the various cell phenotypes coexisting in a tumor is 
known, and the resistance to therapy of each of these regions had 
been estimated, personalized strategies complementary to (or as a 
substitute to) surgery, could be designed to improve chances of 
clinical success. The latter can be understood either as total tumor 
eradication (the standard paradigm as of today) or as achieving 
instead a stable, chronically-controUed tumor burden where less 
aggressive lines keep at bay more resistant ones [102]. In either 
case, significant information towards a personalized treatment 
would be inferred from knowledge of the internal, non-homoge- 
neous structure of a tumor and the resulting differences in therapy 
resistance corresponding to the regions thus identified. 

Unfortunately, to this day only partial information can be 
derived about tumor heterogeneity by means of non-invasive 
medical imaging modalities. Currently available information can 
be mainly used to distinguish various level of oxygen pressure 
within the tumor related to hypoxia processes [9] , necrotic areas or 
highly proliferating regions detected by means of PET techniques 
[103]. While undoubtedly important, such information is not 
enough to design personalized therapies whose results could 
significantly improve those obtained by current procedures. 



In this work, we have proposed a mathematical model of tumor 
growth to gain insight about two key issues: how heterogeneity 
unfolds in a growing tumor, and what type of radiation dosimetry 
is best suited to achieve control in heterogeneous tumors. 
Concerning the first issue, we noticed that substantial information 
about the evolution of spatial heterogeneity within a tumor can be 
retrieved from knowledge of a few key biological properties of the 
tumor cell phenotypes involved. In particular, we have shown in 
the first place that a difference in cell cycle duration between a 
majority of ordinary cancer cells (CCs) and a minority of 
comparatively slow-cycling cancer stem cells (CSCs) leads to a 
concentration of CSCs in regions that can be a priori estimated. In 
the cases just discussed, such regions consist in an internal core 
within an expanding tumor, but the result would apply to other 
geometries as well. In particular, it can be extended to larger 
tumors with corrugated shapes and boundaries. 

We have already mentioned that our key assumption that CSCs 
have longer replication times than CCs is commonplace in the 
literature (see for instance [16], [48], [49], [50], [51], [52], [53], 
[54], [55] and [56]). As a matter of fact, such assumption is 
naturally associated to the consideration of CSCs as a subpopu- 
lation of tumor cells which is able to rescue tumor growth after 
therapies have been delivered. This is related to the fact that 
standard radiation therapies preferentially target dividing cells 
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Figure 8. Time evolution of tumor growth during and after averaged homogeneous radiation therapies. (A) A homogeneous dose of 
2.10 Gy for the case = 0.75 is delivered (Top), and a homogeneous dose of 2.23 Gy for = 0.25 is instead applied (Bottom), assuming in both cases 
of (A) the low migration rate and CSC cycle duration equal to 96 h. (B) A homogeneous dose of 2.36 Gy is delivered (Top) and a homogeneous dose 
of 2.52 Gy (Bottom) for the case /;„ = 0.25 with the high migration rate and CSC cycle durations equal to 96 h (Top) and 48 h (Bottom). In all cases (A, 
B) a standard scheduling (30 sessions along 6 weeks separated by 24 hours intervals except for weekends) was applied. From left to right we show in 
sequential order the tumor before radiotherapy treatment starts, its state after sessions 1 0, 20 and 30, and three stages corresponding to recurrence 
during the period covered (where about 10^ cells is again obtained). Depicted in dark and light green (respectively, dark and light red) are 
proliferating and quiescent CCs (respectively, proliferating and quiescent CSCs). Dead cells are not represented. 
doi:1 0.1 371/journal.pone.0089380.g008 
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Figure 9. Estimates on the total number of CSCs at the end of the recurrence tumor stage. Number of CSCs at the end of the recurrence 
tumor stage (where about 10^ cells is again obtained) and the corresponding standard deviations after performing 20 simulations in each case (with 
different seeds of a random number generator) are shown. (A) For averaged homogeneous therapies corresponding to heterogeneous therapies 
consisting of 2.5 Gy, 2.9 Gy and 3.3 Gy in the inner sphere and 2.0 Gy in the rest of the tumor for the cases y;„ = 0.75, /;„ = 0.50 and = 0.25 (left, 
middle, right) assuming the low migration rate and CSC cycle durations equal to 48 h, 72 h and 96 h (see Table 4). (B) For heterogeneous therapies 
consisting of 2.5 Gy, 2.9 Gy and 3.3 Gy in the inner sphere and 2.0 Gy in the rest of the tumor (Top) and the corresponding averaged homogeneous 
therapies (Bottom) for the cases = 0.75, = 0.50 and ;;„ = 0.25 (left, middle, right) with the high migration rate and CSC cycle durations equal to 
48 h, 72 h and 93 h (see Table 5). In all cases (A, B, C), a standard scheduling (30 sessions along 6 weeks separated by 24 hours intervals except for 
weekends) was applied. Notice that the vertical coordinate is represented in a logarithmic scale. See Tables in the Document SI for further details and 
IVlovies S3, S4, S5, S6, S7 and S8 for some examples of simulations represented in (A), (B) and (C). 
doi:1 0.1 371 /journal.pone.0089380.g009 



(which are more radiosensitive), and thus spare those that have 
slower cycles or remain quiescent. Notice that the cell cycle 
duration could in principle be estimated, at least in vitro, for all cell 
phenotypes known to appear in a given tumor. Importantly, the 
spatial heterogeneity pattern thus observed does not depend so 
much on the precise values of such biological parameters, but 
rather on the fact that they are significantiy diflFerent for the tumor 
cell phenotypes involved. As a consequence, the result obtained is 



robust with respect to fluctuations in cell cycle duration due to 
systemic factors. 

A second result obtained is that, once information about 
functional heterogeneity had been obtained, tumor-tailored 
radiation dosimetries can be designed to improve the treatment 
outcome. We have shown that heterogeneous radiation dosime- 
tries do better than homogeneous ones when regions occupied by 
different radioresistant tumor subpopulations can be identified. 
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Figure 10. Comparing averaged homogeneous radiation therapies without weekend interruptions. Cell survival curves for 20 
simulations (with different seeds of a random number generator) are shown. (A) Averaged dosimetries consisting of AD = 2.W Gy for = 0.75 and 
(B) AD = 2.2'i Gy for /i„ = 0.25, both for the low migration case and CSC cycle duration equal to 96 h. (C, D) Averaged homogeneous therapies 
consisting oi AD = 236 Gy and AD = 2.52 Gy for the case /i„ = 0.25 with the high migration rate and CSC cycle durations equal to 96 h and 48 h 
respectively. The time evolution of CCs and CSCs are represented in green and red respectively. In all cases (A, B, C, D), sessions were scheduled 7 
days a week separated by 24 hours intervals along 30 sessions (without weekend interruptions). Notice that the vertical coordinate is represented in a 
logarithmic scale. See Movie S9 for an example of simulations represented in (D). 
doi:1 0.1 371 /journal.pone.0089380.g010 



and this is in particular the case when more radioresistant 
phenotypes are assumed to replicate at a lower pace than less 
resistant ones. Interestingly, this result holds when the more 
radioresistant phenotype is allowed to sustain unlimited replication 
as in the case of CSCs, as opposed to the limited number of 
replications commonly assumed on CCs. The previous statement 
holds true, no matter the type of scheduling considered (with or 
without weekend interruptions) or the precise result pursued, being 
it total tumor eradication, controlled recurrence or palliative 
treatment. We believe that the comparative advantage of 
heterogeneous radiation dose distributions deserves some consid- 
eration, since to this day homogeneous dosimetries continue to be 
those being commonly implemented worldwide. 

It is worth to stress that our model is quite robust with respect to 
changes in data and parameter values. In particular, our 
conclusions remain in force when CSC and CC cycle durations 
undergo considerable changes, as long as CSC cycle is significantiy 
slower than that of CCs. Also, CCs and CSCs migration rates are 
allowed to undergo substantial changes (corresponding for 
instance to slow and fast migration processes) as far as both 
tumor cell phenotypes share a similar migration rate. Moreover, 
our results continue to hold when changes in the choice of the 
radiosensitivity parameters a. and /? in (1) are allowed, or when 
different fractions of the minority phenotype (CSC) are assumed. 
For instance, our results are not confined to the choice made for 



the assumed percentage (15%) of CSCs present at an early stage of 
tumor growth. They continue to hold if a different figure for that 
proportion is taken, as long as CSCs remain a small fraction of the 
total tumor population at that stage. 

On the other hand, cancer cell plasticity has recentiy received 
considerable attention ([97], [98] and [99]). We have studied a 
particular example where in addition to CCs being generated by 
CSCs with asymmetric division, a small percentage of CCs are 
transformed to CSCs as a consequence of radioresistance to 
therapy. Our conclusion that heterogeneous, tumor-adapted 
radiation therapies fare better than their corresponding averaged 
homogeneous versions continues to hold in this case. 

We conclude by discussing on some of the limitations of this 
work, as well as on possible extensions thereof. To begin with, we 
are aware that more research is needed to understand the possible 
mechanisms that can be responsible for slow cycling of CSCs. 
Particularly relevant in this context would be to ascertain if slow 
cycling can, at least in some cases, be established as an intrinsic 
property of CSCs or if it could alternatively be induced by systemic 
feedback in the course of tumor growth. Interestingly, even if 
CSCs are assumed to cycle faster than CCs, our model still shows 
that heterogeneous dosimetries adapted to the resulting tumor 
heterogeneity continue to outperform standard homogenous 
therapies currently in use. 
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Figure 11. Comparing the effects of lower radiation dosimetries with and without hyperfractionation. Cell survival curves for 20 
simulations (with different seeds of a random number generator) are shown in the case/;,, = 0.25 for the low migration case and CSC cycle duration 
equal to 96 h. (Top) From left to right heterogeneous therapies consisting of 2.3 Gy (A) and 1.7 Gy (D) in the inner sphere, and 1.8 Gy (A) and 1.2 Gy 
(D) in the rest of the tumor respectively. (Middle) From left to right the averaged homogeneous therapies corresponding to 1 .9 Gy (B) and 1 .3 Gy (E) 
are represented. Radiation dose delivery been made according to 5 days a week, 30 sessions in total, at 24 hours (A, B) and at 12 hours (D, E) intervals 
with weekend interruptions. The time evolution of CCs and CSCs is represented in green and red respectively. (Bottom) Number of CSCs and the 
corresponding standard deviations at the end of the recurrence tumor stage (where about 1 0^ cells is again obtained) for heterogeneous (yellow) and 
averaged homogeneous (blue) radiation therapies (C, F). Notice that the vertical coordinate is represented in a logarithmic scale. See IVIovie SI 0 for an 
example of simulations represented in (E). 
doi:1 0.1 371 /journal.pone.0089380.g01 1 



A general conclusion that follows from our study is that detailed 
information about intratumoral heterogeneity is needed in order 
to implement efficient dose-painting techniques in clinical practice. 
In particular, in this work a clear dependence on tumor 
heterogeneity of the radiation doses needed to achieve tumor 
control has been obtained. In fact, the inner tumor regions where 
more radioresistant tumor cell phenotype remains confined are 
shown to strongly depend on CSC cycle duration and their 
probability of asymmetric division. In the particularly unfavorable 



assumption of fast CSCs cycling, this region may rank from 20% 
to 100% of the total tumor volume. In this latter situation, a worst- 
case scenario corresponding to a high and homogeneous radiation 
dose being prescribed and only limited by neighboring organs at 
risk tolerance, is recovered that corresponds to current clinical 
practice. Our results suggest that such situation could be 
considerably improved in many cases if and when sufficient 
information about key different biological and radiobiological 
properties of the tumor cell phenotypes present in a given tumor is 
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available, be it either by estimating patient-specific parameters or 
by means of medical imaging techniques. 

Finally, it looks feasible from a mathematical viewpoint to 
address within this framework situations where larger tumors are 
considered, a number of cell phenotypes coexist there due to 
mutations, and other effects (immune response, nutrient limitation, 
etc.) are accounted for. For example, the modeling framework 
selected in this work permits simulations to be scaled up to cubic 
centimeter sizes, though at the expense of lower spatial and 
functional resolution, and more computing resources. It can also 
be used to construct hybrid models, zooming in at the cell scale in 
regions of interest. In particular, we have chosen to represent each 
cell individually to exclude averaging effects when studying the 
relation between tumor heterogeneity and simulated radiation 
outcomes. We hope that this work could provide a starting point 
towards the study of the more general situations described above. 

Supporting Information 

Document SI Details of computer simulations of the 
model of tumor growth and additional results. 

(PDF) 

Movie SI Time evolution of a homogeneous tumor 
growth where only CCs are present for the low 
migration case. Tumor growth is allowed unchecked from a 
size of about 10^ cells until about 10^ cells are present. Then, a 
standard homogeneous radiation therapy corresponding to 30 
sessions of 20 Gy each is delivered. Sessions from Monday to 
Friday are scheduled along 6 weeks, separated by 24 hours 
intervals except for weekends, where a 72 hours inter\'al is allowed. 
Depicted in dark and light green (respectively, dark and Ught red) 
are proliferating and quiescent CCs (respectively, proliferating and 
quiescent CSCs). Dead cells are represented in black. See Figure 3 
(A, B) in the article. 
(MP4) 

Movie S2 Time evolution of a homogeneous tumor 
growth where only CCs are present for the high 
migration case. Tumor growth is allowed unchecked from a 
size of about 10 ' cells until about 10*' cells are present. Then, a 
standard homogeneous radiation therapy corresponding to 30 
sessions of 2.0 Gy each is delivered. Sessions from Monday to 
Friday are scheduled along 6 weeks, separated by 24 hours 
intervals except for weekends, where a 72 hours inter\'al is allowed. 
Depicted in dark and light green (respectively, dark and light red) 
are proliferating and quiescent CCs (respectively, proliferating and 
quiescent CSCs). Dead cells are represented in black. See Figure 3 
(C, D) in the article. 
(MP4) 

Movie S3 Time evolution of a heterogeneous tumor 
grovtrth (where CCs and CSCs are present) for the low 
migration case with Pa = 0J5 and CSC cycle duration 
equal to 96 h. Tumor growth is allowed unchecked from a size 
of about 10'^ cells until about 10** cells are present. Then, a 
homogeneous radiation therapy consisting of 2.10 in the tumor 
is delivered. Treatment sessions are scheduled along 6 weeks, 
separated by 24 hours intervals except for weekends, where a 72 
hours interval is allowed. Depicted in dark and light green 
(respectively, dark and light red) are proliferating and quiescent 
CCs (respectively, proliferating and quiescent CSCs). Dead cells 
are r(;presented in black. See Figure 7 (B) in the article. 
(MP4) 

Movie S4 Time evolution of a heterogeneous tumor 
grovtrth (where CCs and CSCs are present) for the low 



migration case with /Jo = 0.25 and CSC cycle duration 

equal to 63 h. Tumor growth is allowed unchecked from a size 
of about 10^ cells until about 10® cells are present. Then, a 

homogeneous radiation therapy consisting of 2.23 Gji in the tumor 
is delivered. Treatment sessions are scheduled along 6 weeks, 
separated by 24 hours int(-r\ als cxcc'iit for weekends, where a 72 
hours interval is allowed. Depicted in dark and light green 
(respectively, dark and light red) are proliferating and quiescent 
CCs (respectively, proliferating and quiescent CSCs). Dead cells 
are represented in black. See Figure 7 (D) in the article. 
(MP4) 

Movie S5 Time evolution of a heterogeneous tumor 
growth (where CCs and CSCs are present) for the low 
migration case with pg = 0.25 and CSC cycle duration 

equal to 48 h. Tumor growth is allowed unchecked from a size 
of about 10^ cells until about 10® cells are present. Then, a 
heterogeneous radiation therapy consisting of 3.3 Gf in the inner 
sphere and 2.0 Gy in the rest of the tumor is delivered. Treatment 
sessions are scheduled along 6 weeks, separated by 24 hours 
intervals except for weekends, where a 72 hours interx^al is allowed. 
Depicted in dark and light green (respectively, dark and light red) 
are proliferating and quiescent CCs (respectively, proliferating and 
quiescent CSCs). Dead cells are represented in black. See Tables 2 
and 4 in the article for further details. 
(MP4) 

Movie S6 Time evolution of a heterogeneous tumor 
growth (where CCs and CSCs are present) for the low 
migration case with Pa = 0.25 and CSC cycle duration 
equal to 48 h. Tumor growth is allowed unchecked from a size 
of about 10^ cells until about 10® cells are present. Then, a 
homogeneous radiation therapy consisting of 2.32 Gj) in the tumor 
is delivered. Treatment sessions are scheduled along 6 weeks, 
separated by 24 hours intervals except for weekends, where a 72 
hours inter\'al is allowxxi. D("picted in dark and light green 
(respectively, dark and light red) are proliferating and quiescent 
CCs (respectively, proliferating and quiescent CSCs). Dead cells 
are represented in black. See Tables 2 and 4 in the article for 
further details. 
(MP4) 

Movie S7 Time evolution of a heterogeneous tumor 
growth (where CCs and CSCs are present) for the high 
migration case with pa = 0.25 and CSC cycle duration 

equal to 48 h. Tumor growth is allowed unchecked from a size 
of about 10^ cells until about 10® cells are present. Then, a 
heterogeneous radiation therapy consisting of 3.3 Gy in the inner 
sphere and 2.0 Gf in the rest of the tumor is delivered. Treatment 

sessions are scheduled along 6 weeks, separated by 24 hours 
intervals except for weekends, where a 72 hours interval is allowed. 
Depicted in dark and light green (respectively, dark and light red) 
are proliferating and quiescent CCs (respectively, proliferating and 
quiescent CSCs). Dead cells are represented in black. See Figure 7 
(G) in the article. 
(MP4) 

Movie S8 Time evolution of a heterogeneous tumor 
growth (where CCs and CSCs are present) for the high 
migration case with Pa = 0.25 and CSC cycle duration 
equal to 48 h. Tumor growth is allowed unchecked from a size 
of about 10^ cells until about 10® cells are present. Then, a 
homogeneous radiation therapy consisting of 2.52 Oyin the tumor 
is delivered. Treatment sessions are scheduled along 6 weeks, 
separated by 24 hours intervals except for weekends, where a 72 
hours interval is allowed. Depicted in dark and light green 
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(respectively, dark and light red) are proliferating and quiescxnt 
CCs (respectively, proliferating and quiescent CSCs). Dead cells 
are represented in black. See Figure 7 (H) in the article. 
(MP4) 

Movie S9 Time evolution of a heterogeneous tumor 
growth (where CCs and CSCs are present) for the high 
migration case with />a = 0.25 and CSC cycle duration 

equal to 48 h. Tumor growth is allowed unchecked from a size 
of about 10^ cells until about 10^ cells are present. Then, a 
homogeneous radiation therapy consisting of 2.52 Gy in the tumor 
is delivered. Treatment sessions are scheduled 7 days a week 
separated by 24 hours intervals along 30 sessions (without weekend 
interruptions). Depicted in dark and light green (respectively, dark 
and light red) are proliferating and quiescent CCs (respectively, 
proliferating and quiescent CSCs). Dead cells are represented in 
black. See Figure 10 (D) in the article. 
(MP4) 

Movie SIO Time evolution of a heterogeneous tumor 
grovtrth (where CCs and CSCs are present) for the low 
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